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1.0  Introduction 
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Motivation 


There  has  been  a  considerable  amount  of  work  done  on  understanding,  modeling  and  computing 
the  observed  self  focusing  and  ionization  of  ultra  short  laser  pulses  (USP)  when  they  propagate 
in  air  and  other  materials.  There  is  a  need  for  study  of  the  potential  applications  of  these  effects. 
One  application  that  has  been  suggested  is  the  use  of  a  USP  ionization  channel  to  guide  and 
extend  the  length  of  high  voltage  electrical  arcs.  Although  there  are  questions  about  the  feasible 
utility  of  this  idea  (Ref  1),  we  need  to  develop  new  physics  models  and  analysis  techniques  for 
assessing  this  utilization. 

Background  Discussion 

Modeling  of  USP  laser  pulse  propagation  with  self  focusing  and  ionization  has  been  done 
previously  Braun  (Ref  2),  Mlejnick  (Ref  3),  and  Sprangle  (Ref  4).  A  study  was  done  at 
AFRL/DE  in  order  to  understand  the  reported  associated  EMP  radiation.  Analytical  estimates 
were  discussed  and  a  numerical  atmospheric  propagation  code  was  developed,  based  on  standard 
pulse  envelope  techniques  outlined  in  (Ref  4).  The  modeling  effort  also  included  development  of 
a  full  Maxwell’s  Equation  model  for  propagation  and  material  interaction.  This  modeling  and 
computational  simulation  development  has  been  published  in  (Ref  1). 

The  possibility  of  electrical  energy  discharge  guided  by  an  USP  produced  ionization  channel  in 
air  has  been  raised.  As  an  example  the  positive  tenninal  of  a  high  voltage  source  may  be 
connected  to  a  conductor  in  close  proximity  to  the  filament  produced  in  air  by  an  ultra  short  laser 
and  the  return  electrode  may  consist  of  the  target,  and  the  material  that  the  target  is  standing  on. 
Although  the  use  of  Ultra-Violet  (UV)  lasers  to  trigger  electrical  breakdown  is  well-known,  ultra 
short  laser  produced  filaments  would  appear  to  offer  an  advantage  in  precise  control  of  the 
location  and  direction  of  the  electrical  breakdown.  Experiments  performed  to  date  have  shown 
some  promise  but  a  satisfactory  modeling  capability  is  necessary  to  guide  experiments  and  to 
assess  the  viability  of  potential  practical  uses. 

The  scope  of  this  project  included  the  formulation  of  an  overall  simulation  methodology  and  the 
development,  testing  and  integration  of  the  requisite  sub  models.  It  did  did  not  include  porting  to, 
and  parallel  operation  on,  a  High  Performance  Computer  (HPC)  platform.  This  will  be  required 
to  determine  the  ultimate  length  of  streamer  propagation  and  the  energy  delivery  to  an  anode 
for  cases  of  practical  application. 

Problem  Overview 


In  the  current  effort  an  electrical  discharge  model  was  developed  to  describe  the  USP  induced 
breakdown  of  a  spark  gap  at  ambient  conditions.  The  problem  geometry  is  shown  in  Figure  1. 
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The  USP  induced  ionization  is  schematically  represented  by  the  arrows  between  the  two  disk 
electrodes. 


Figure  1.  Schematic  view  of  the  electrical  breakdown  problem. 


The  problem  of  predicting  the  time  behavior  of  an  electrical  discharge  has  been  widely  studied 
by  many  investigators  and  is  still  not  fully  understood.  Examples  of  electrical  discharges  range 
from  the  local  corona  discharge  induced  on  high  voltage  transmission  lines,  to  the  long  spark 
spanning  a  gap  with  dimensions  of  centimeters  to  tens  of  meters,  to  the  phenomena  of  natural 
lightning.  The  problem  is  recognized  to  be  difficult  due  to  the  large  range  of  dimensions  and 
time  scales  involved,  as  well  as  an  incomplete  understanding  of  all  the  relevant  physics. 

The  physics  of  electrical  breakdown  in  a  one-dimensional  planar  geometry  was  studied  by  Kline 
[Ref.  5]  using  Monte-Carlo  techniques.  The  velocities  of  anode  and  cathode  directed  streamers 
in  nitrogen  at  relatively  low  pressures  were  compared  to  experiment.  Of  particular  interest  was 
the  use  of  the  empirical  results  of  Penny  [Ref.  6]  to  determine  the  photo  ionization  rates  in 
nitrogen.  Monte  Carlo  methods  have  also  been  demonstrated  by  Kunhardt  [Ref.  7]  and  others. 
The  Monte  Carlo  approach  has  the  obvious  advantage  of  being  able  to  represent  the  six 
dimensional  distribution  function  of  the  electrons  and  ions.  Flowever  the  method  is 
computationally  intense  as  the  random  velocity  of  the  electron  distribution  has  to  be  properly 
advanced  over  a  spatial  grid  with  relatively  small  length  scales.  As  a  result  much  effort  has  been 
put  into  developing  models  that  effectively  integrate  over  the  distribution  function  and  advance 
the  electron  and  ion  drift  velocities.  These  approaches  either  use  the  electron  and  ion  mobilities 
by  Yoshida  [Ref.  8],  Wu  [Ref.  9],  Wang  [Ref.  10],  Vitello  [Ref.  1 1],  Kulikovsky  [Ref.  12],  or 
solve  the  electron  swami  equations  by  Guo  [Ref.  13]. 

In  the  calculations  described  here  the  ionization  path  produced  by  the  self  focusing 
propagation  of  ultra  short  laser  pulses  was  included  as  a  fixed  (time  independent  pre-ionization). 
Time  dependant  laser  ionization  during  the  formation  of  the  discharge  channel  can  easily  be 
added.  The  levels  of  the  USP  laser  ionization  used  were  taken  from  previous  studies  of  self 
focusing  propagation  by  Zimmerman,  et  al  [Ref.  1]. 
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The  use  of  electron  and  ion  mobilities  to  relate  the  electron  and  ion  velocities  directly  to  the  local 
electric  field  strength  is  commonly  referred  to  as  the  equilibrium  approach  as  it  assumes  the 
electron  and  ion  have  come  into  equilibrium  between  the  effective  collision  frequencies  and  the 
acceleration  of  the  electric  field.  The  use  of  a  full  set  of  fluid  equations  to  model  the  electrons 
and  ions  is  referred  to  as  the  non-equilibrium  approach  and  is  more  demanding  on  computer 
resources  and  may  be  more  susceptible  to  instabilities  when  the  electron  plasma  frequency 
approaches  the  inverse  of  the  time  step.  Although  neither  approach  has  the  fidelity  of  a  full 
Monte  Carlo  simulation  of  the  electron  and  ion  distribution  functions  both  the  equilibrium  and 
non-equilibrium  approaches  are  felt  to  be  adequate  for  electrical  breakdown  problems  due  to  the 
simple  fact  that  the  electrons  and  ions  travel  “one  way”  between  the  cathode  and  anode  (or  visa- 
versa). 

The  majority  of  papers  referenced  above  consider  breakdown  in  nitrogen  and  do  not  consider  late 
time  heating  and  arcing  effects.  Arcing  by  itself  has  been  studied  by  Plooster  [Ref.  14]  and 
Ganesh  [Ref.  15]  but  the  models  used  are  typically  one  dimensional  radial  models  and  consider 
thennal  and  hydrodynamic  effects  at  relatively  late  times  as  compared  to  the  streamer  breakdown 
problem.  There  has  been  more  recent  interest  in  modeling  electrical  breakdown  in  air  and 
incorporating  a  more  complete  description  of  the  complex  air  chemistry;  see  Kossyi  [Ref.  16], 
Naidis  [Ref.  17],  Aleksandrov  [ref.  18].  Many  of  the  simulations  are  PA  D  where  the  streamer 
radius  is  fixed.  The  behavior  of  a  long  leader  in  air  has  been  studied  using  ID  radial  models 
Aleksandrov,  [Ref.  19]  and  [Ref.  207]  and  a  PA  D  model,  Aleksandrov  [Ref.  21].  The  PA  D  or 
quasi  2D  simulations  are  often  preferred  as  the  spatial  grid  requirements  and  the  need  to  solve 
the  stiff  system  of  complex  air  chemistry  makes  the  problem  computationally  challenging.  A 
contemporary  review  of  the  subjects  of  spark  discharge  and  air  chemistry  can  be  found  in 
Bazelyan  and  Raizer  [Ref.  22]  and  Capitelli  [Ref.23]. 

In  this  report  we  will  be  interested  in  the  problem  of  the  initiation  of  an  electrical  discharge  by 
means  of  a  laser  induced  ionized  path  between  two  oppositely  charged  electrodes  in  air  at 
atmospheric  pressure.  In  addition  we  are  interested  in  both  the  early  time  streamer  dynamics  as 
well  as  the  development  of  a  hot  electrically  conductive  leader  channel  that  allows  the  electrical 
breakdown  to  travel  over  large  distances.  The  problem,  as  illustrated  in  Figure  1,  consists  of  a 
spark  gap  constructed  of  two  cylindrical  electrodes  separated  by  a  distance  d  .  The  spark  gap  is 
driven  by  a  voltage  source  Vs  with  a  series  impedance  Rs  resulting  in  voltages  Vu ,  V ,  on  the  upper 
and  lower  plates.  A  laser  is  used  to  create  a  channel  of  ionization  along  the  axis.  The  problem  is 
to  predict  the  time  history  of  the  electrical  discharge  including  the  discharge  current  I  and  the 

voltage  drop  Vg=Vu- V,  across  the  electrodes. 

Modeling  Approach 

The  modeling  approach  described  here  divides  the  electrical  breakdown  problem  into  a 
sequential  set  of  calculations  based  on  the  time  scale  and  the  physical  mechanisms  of  the 
problem.  Specifically  we  will  divide  the  problem  into  an  early  time  and  a  late  time  regime. 


Early  Time  Phase:  The  initial  electrical  breakdown  problem  consists  of  the  electron  avalanche 
resulting  in  space  charge  waves.  Depending  on  the  details  of  the  electrodes  and  the  use  of  a  laser 
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for  pre-ionization  it  is  necessary  to  consider  electron  transport,  ion  transport  and  diffusion.  In 
this  report  the  early  time  electron  motion  will  be  modeled  with  a  fluid  approach.  It  will  be 
necessary  to  consider  various  source  terms  including  avalanche,  electron  attachment  and 
electron-ion  recombination.  During  the  early  time  phase  of  the  electrical  discharge  the  electron 
temperature  will  not  be  in  equilibrium  with  the  ion  and  neutral  temperatures. 


Late  Time  Phase:  During  the  later  part  of  the  breakdown  process  the  neutral  gas  is  expected  to 
become  hot  enough  so  that  thermal  agitation  accounts  for  most  of  the  free  charge.  In  this  case 
we  will  assume  thermal  equilibrium  between  the  electrons,  ions  and  neutral  gas  molecules  and 
atoms  and  use  an  electrical  conductivity  to  model  the  final  breakdown  process. 

A  two  dimensional  numerical  model  in  cylindrical  coordinates  will  be  developed  to  model  the 
electrical  discharge  problem.  In  addition  a  one-dimensional  model  will  be  developed.  The  one¬ 
dimensional  model  is  relatively  easy  to  understand  and  so  is  useful  for  developing  and  testing 
numerical  methods  and  physical  assumptions. 

Simulating  a  Laser  Induced  Electrical  Discharge 

The  problem  of  interest  is  determining  whether  a  laser  trigger  electrical  discharge  can  propagate 
over  a  distance  of  several  meters.  As  will  be  shown,  the  early  time  breakdown  problem  will 
require  -5-10  micron  resolution  over  the  dimensions  of  the  head  of  the  discharge.  This 
dimension  is  on  the  order  of  centimeters.  As  a  result  it  is  not  possible,  nor  necessary  to  actually 
simulate  the  entire  discharge  process.  Instead  the  approach  taken  in  simulating  the  nuclear 
induced  lightning  problem  by  Gardner  [Ref.24]  will  be  followed.  To  recall  the  nuclear  lightning 
discharge  involves  the  propagation  of  a  streamer  in  the  electrically  conductive  atmosphere 
induced  by  the  prompted  and  delayed  gamma  radiation  from  a  nuclear  detonation.  The  approach 
used  was  to  calculate  the  local  discharge  physics  in  a  region  of  about  10  cm  at  the  tip  of  a 
conductive  streamer.  If  the  numerical  model  showed  that  the  discharge  would  propagate,  rather 
than  die  out,  then  the  simulation  was  successful.  If  necessary  the  simulation  may  be  repeated  for 
various  streamer  lengths. 

Model  Geometries 


A  two  dimensional  simulation  model  was  developed  to  model  the  propagation  of  a  laser  guided 
arc  over  short,  centimeter  scale,  ranges.  A  secondary  one  dimensional  model  was  also  used  to 
facilitate  development  and  implementation  of  various  physics  sub  models  required  for  the 
complete  simulation.  The  range  limitation  results  from  run  time  feasibility  on  single  PC 
workstation  computer.  This  range  can  be  extended  to  the  10’s  of  cm  region  by  porting  the  model 
to  a  HPC.  This  procedure  has  been  investigated  and  outlined  but  could  not  be  done  within  the 
resource  limitations  of  this  effort.  A  further  method  was  researched  that  will  allow  the 
simulation  to  be  used  over  many  meter  ranges  by  utilizing  a  stepped  set  of  local  -10  cm 
snapshot  calculations  at  different  ranges;  these  can  be  used  to  approximate  propagation  over  the 
full  streamer  length.  The  simulation  is  based  on  a  cylindrically  symmetric  computational  volume 
that  contains  metallic  objects  (boundary  conditions).  The  simulation  is  used  to  predict  the 
breakdown  of  the  spark  gap  electrode  geometry  shown  in  Figure  2.  The  spatial  grid  for  the 
ionization  is  limited  to  points  near  the  center  ionized  column  while  the  electrostatic  grid  covers 
the  entire  problem  space. 
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Figure  2.  Simulation  geometry  used  for  numerical  simulation  of  the  electrical 
breakdown  and  ionization  guided  arc. 


Computational  Overview 

The  ionized  air  is  represented  as  a  5  species  model 


-^-  =  -{C  +  B)n+  Enm  +  Rjint  +  Run+n_  - Se 

+  V  {neUe  )  -  DeV\  =  (C  "  A)  "e  ~  ReW+  +  Se 
_ 

~^T' +  V '  ( n+u+  )  =  Cne  -  Rein+ne  -  Run+n_  +  Se  ( 1 ) 

^  _ 

— -  +  V  •  in  v  )  =  Ana  -  Rjiji 

dt  y  ~  ~ 


dt 


=  Bne 


-En 


m 


The  variables n0,ne,n+,n_,nm  refer  to  the  neutral  air,  electron,  positive  ion,  negative  ion,  and 

excited  neutral  molecule  densities.  The  source  terms  include  the  Ultra-Violet  (UV)  ionization  of 
the  air  and  the  cathode  via  the  inelastic  neutral  molecule  densities.  These  five  equations  are 
simultaneously  solved  in  each  finite  difference  cell  using  implicit  integration  techniques.  The 
electric  field  is  determined  using  a  two  dimensional  FFT  technique  [Ref.  22,  23]  and  the 
electron/ion  velocities  are  advanced  using  mobilities.  The  entire  system  of  equations  is  then 
iterated  with  the  circuit  equations  to  achieve  a  self  consistent  solution  with  the  voltage  source 
and  load  impedance. 
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Arc  formation  requires  the  heating  of  the  air.  The  macroscopic  conductivity  of  the  heated  air  can 
be  use  to  replace  the  detailed  microscopic  treatment  of  charge  flow  at  later  times  in  the 
simulation.  The  change  in  specific  energy  of  neutral  air  is 


d[ps) 

dt 


=  E  -J  +  <JepE 


4  °hr 

iP 


(2) 


The  temperature  and  electrical  conductivity  were  then  determined  from  Plooster’s  equation  of 
state  and  electron  mobility  model  [Ref.  15],  [Ref.  27]  given  by 


e(P’T)~  ~(5  +  4>)+3(4 +40 


RT  ~\~  A  I  +  A,l \  +  Arylry 

o  o  11  2  2 


(3) 


where  e  is  the  specific  internal  energy,  A0  represents  the  fraction  of  the  air  molecules  that  are 
dissociated  and  Al,A2  represent  the  fraction  of  air  atoms  that  have  been  signally  or  doubly 
ionized.  These  coefficients  are  functions  of  temperature  T  and  density  p  .  The  electron  density 
due  to  thennal  ionization  is  given  by 

neP=2PrN0(A+  Al)  .  (4) 

From  the  change  in  internal  specific  energy,  it  is  possible  to  numerically  solve  the  equation  of 
state  for  temperature.  The  electron  mobility  is  given  by 


Pen{m2lV~S) 


2.64 

Pr(l-Ai)jT{K) 


2.45x1  O'2 


(5) 


The  use  of  Plooster’s  equation  of  state  eliminates  the  requirement  to  solve  a  potentially  very 
large  (-100)  system  of  air  chemistry  equations. 


Poisson  Computation  Timing 


The  electric  field  solution  time  is  dominated  by  the  solution  of  Poisson’s  equation  that  finds  a 
quasi-static  field  solution  due  to  a  distribution  of  charged  particles.  The  fast  Poisson  solver  that 
we  have  developed  fits  the  grid  with  splines  in  one  (radial)  dimension  and  transforms  the  other 
free  dimension  (z)  via  a  sine  transform  to  rapidly  solve  the  differential  equation  in  the 
“frequency  domain”  followed  by  an  inverse  transform  back  to  the  original  spatial  domain.  The 
dominate  sub-computation  in  the  solution  of  the  electric  field  is  the  forward  and  inverse  sine 
transfonns.  The  procedure  for  parallelizing  FFTs  is  well-documented  [Ref.  28]. 


Figure  3  shows  computational  timing  results  comparing  the  FFT-based  fast  Poisson  solver, 
developed  for  this  simulation,  to  a  Successive  Over-Relaxation  (SOR)  method.  For  this  test 
case,  the  potential  is  solved  for  a  conducting  sphere  at  the  origin  of  a  square  computational  grid. 
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The  red  curve  shows  the  computation  time  for  the  fast  solver;  the  blue  the  computation  time  for 
the  SOR  method;  each  as  a  function  of  grid  size.  The  pink  and  green  lines  are  present  merely  to 
guide  the  eye:  they  approximately  match  the  curves  at  the  first  point  and  show  the  expected  times 
for  algorithms  that  scale  as  the  square  of  the  linear  dimension.  That  is,  they  show  an  ideal 
scaling  for  a  problem  that  grows  in  proportion  to  the  number  of  cells.  The  fast  FFT  method  does 
a  reasonable  job  at  following  the  curve  up  until  the  last  point  (8196  points)  where,  most  likely, 
the  size  of  the  arrays  exceed  the  cache  of  the  processor  and  an  additional  time  is  required  due  to 
swap  memory  in  and  out  of  the  CPU.  If  the  supposition  is  true,  the  problem  should  be  negated 
by  our  proposed  parallel  approach.  The  time  of  the  SOR  computation  grows  so  quickly  that  it 
did  not  make  sense  to  compute  a  solution  for  a  grid  larger  than  1024  by  1024. 


Comparison  between  FFT  and  SOR  Poisson  Solvers 


Figure  3.  Timing  Comparison  of  FFT  and  SOR  Based  Methods. 


Considerations  for  Utilizing  Parallel  Processing 

Though  beyond  the  resources  of  the  current  effort  the  mechanics  of  performing  the  conversion  to 
parallel  processing  have  been  studied  and  outlined.  This  extension  can  be  accomplished  using  the 
standard  Message  Passing  Interface  (MPI.)  In  particular,  the  Local  Area  Multicomputer- 
Message  Passing  Interface  (LAM/MPI)  can  use  the  MPI  implementation  [Ref.  29].  MPI 
provides  Fortran  and  C  routines  that  allow  easy  and  safe  communication  between  processors.  To 
model  electrical  breakdown  for  geometries  where  the  electrons  are  separated  by  a  meter  or  more 
it  will  be  necessary  to  break  the  problem  up  into  a  series  of  smaller  snapshots.  Each  snapshot 
will  model  the  streamer  head  over  dimensions  of  approximately  10  cm  at  various  locations 
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between  the  cathode  and  the  anode.  If  the  discharge  is  seen  to  propagate  for  each  shapshot  then 
clearly  the  discharge  will  propagate  from  the  cathode  to  the  anode.  A  similar  approach  has  been 
used  for  the  Nuclear  Lightning  problem  [Ref.21].  This  process  is  illustrated  in  Figure  4.  The  left 
hand  part  of  the  Figure  shows  the  relationship  between  the  larger  finite  difference  spatial  grid, 
where  the  figure  shows  the  relationship  between  the  larger  finite  difference  spatial  grid,  where 
the  electrostatic  solution  is  solved,  and  one  particular  location  of  the  sub-grid  where  the  rate 
equations  are  advanced.  The  USP  induced  ionization  channel  is  located  along  the  axis  between 
the  cathode  and  the  anode  and  the  discharge  moves  from  bottom  to  top  in  Figure  4.  below. 

The  sub-grid,  the  electrical  discharge  is  known  to  be  a  hot  leader  below.  As  such  it  must  have  a 
relatively  high  electrical  conductivity  and  its  potential  is,  therefore,  set  to  the  same  potential  as 
the  cathode.  The  right  hand  part  of  the  Figure  shows  the  propagation  of  the  electrical  discharge 
within  the  sub-grid.  The  decision  as  to  whether  the  discharge  is  propagating  through  the  sub-gird 
is  straightforward.  If  the  shape  and  velocity  of  the  tip  of  the  discharge  is  the  same  at  two  times 
along  the  path,  say  Tstart  and  Tstop,  then  the  discharge  is  propagating  and  not  dissipating.  In 
that  case,  the  projected  discharge  path  over  the  entire  extent  of  the  sub-grid  is  set  to  the  electrical 
conductivity  of  the  hot  leader  and  the  location  of  the  rate  equation  sub-grid  is  moved  up  by 
approximately  10  cm  along  the  discharge  path  as  the  leader  propagates.  This  process  may  be 
repeated  over  the  entire  cathode  anode  separation  path  so  that  the  sub-grid  is  co-located  with  the 
tip  of  the  leader. 


Tstop 


Tstart 


Figure  4.  Scheme  for  implementing  a  snapshot  calculation  for  extended  ranges. 
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2.0  Theoretical  Formulation 


The  fluid  approach  for  electrons  can  be  carried  out  by  equilibrium  or  a  non-equilibrium  model. 

In  both  cases  the  quantities  of  interest  are  the  number  densities  and  the  drift  velocities.  In  the 
equilibrium  model  it  is  assumed  that  the  acceleration  due  to  the  local  electric  field  is  balanced  by 
the  collision  rate  and  so  the  drift  velocity  is  a  function  of  the  local  electric  field.  However  it  is 
not  assumed  that  the  electrons  and  ions/neutrals  have  the  same  temperature  and  during  the  initial 
breakdown-streamer  phase  we  expect  that  the  electron  temperature  is  significantly  higher  than 
the  ion/neutral  temperature  i.e.  Te  »T.,T0.  In  the  non-equilibrium  model  the  drift  velocity  is 
determined  by  solving  the  fluid  equation  for  momentum.  In  the  equilibrium  model  we  need  the 
continuity  equation  for  the  electrons  and  ion  species,  electric  field  mobilities  to  relate  the 
electron  and  ion  velocities  to  the  local  electric  field,  and  electric  field  dependent  rate  coefficients 
for  the  various  processes  that  add  and  subtract  to  the  various  number  densities.  In  the  non¬ 
equilibrium  model  we  need  continuity  equations  for  number  density,  momentum,  and  energy 
density,  as  well  as  an  equation  of  state.  In  this  case  the  rate  coefficients  are  typically  written  in 
tenns  of  energy. 

Although  the  non-equilibrium  approach  offers  a  more  detailed  description  of  the  breakdown 
dynamics  the  solution  of  the  fluid  equations  adds  considerable  complexity  to  the  problem. 
Therefore  the  equilibrium  approach  will  be  initially  followed.  The  fluid  equations  needed  for  the 
non-equilibrium  approach  are  separately  documented. 


Equilibrium  Model 

The  equilibrium  model  for  the  electron  and  ion  species  is  given  by  a  set  of  continuity  equations 
i.e. 


=  Sources  -  Sinks 
dt 

Qyi  _ 

— -  +  V  •  ( n  v„ )  -  Dy2n ,  =  Sources  -  Sinks 
8t  v  ’ 

^  _ 

— —  +  V  •  ( n  U. )  =  Sources  -  Sinks 
dt  K  +  +) 

_ 

— -  +  V  •  (n  v  )  =  Sources  -  Sinks 
dt  v 

.  =  Sources  -  Sinks. 


(6) 


dt 


plus  a  set  of  mobility  equations 
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(V) 


Ve=~Me(\E\)E 
v+=np^E^E 
o_=-Vp{\E\)E. 

In  the  continuity  equations  nQ  is  the  number  density  of  neutral  molecules  in  the  ground  state,  ne  is 
the  number  density  of  electrons,  n+  is  the  number  density  of  singly  charged  ions,  n  is  the 
number  of  negatively  charged  ions,  and  nm  is  the  number  density  of  neutral  molecules  in  an 
average  excited  state.  In  general  the  source  and  sink  terms  can  have  an  arbitrarily  large  number 
of  cross  coupling  tenns.  Note  that  we  are  dropping  the  divergence  term  for  the  neutral  molecules 
and  the  diffusion  term  for  all  the  ion  species. 


The  various  processes  contained  in  the  source  and  sink  terms  include  electron  avalanche, 
electron  attachment,  electron-ion  and  ion-ion  recombination,  inelastic  excitation  and  decay  of 
neutral  molecules,  and  photo-ionization  of  neutral  molecules.  As  a  result  we  can  rewrite  the 
continuity  equations 


8n„ 

dt 

dne 

dt 

dn+ 

dt 

dn_ 

dt 

dn 

m 

dt 


=  -(C  +  B)ne+  Enm  +  Reinen+  +  Run+n_  -  Se 
+  V  •  (neDe )  -  DeV2ne  =(C-A)ne-  Reinen+  +  Se 
+  V-(n+U+)  =  Cne  - Rein+ne  - Run+n  +  Se 
+  V  •  [n  o  )  =  Ane  -  Run+n_ 

=  Bne  —Enm. 


(8) 


These  equations  can  be  written  in  the  general  fonn 


dnj 

dt 


X a,  ox +XXa  O'.  k)njnk  + 

k  j  k 


(9) 


We  can  define  an  early  and  a  late  time  regime  for  the  calculations.  The  early  time  regime 
models  the  propagation  of  the  avalanche  and  streamer.  At  late  times  when  the  electron 
avalanche  has  crossed  the  gap  between  the  two  electrodes  the  problem  should  settle  to  a  steady 
state  where  the  gradients  of  the  electron  and  ion  densities  become  relatively  constant  in  time.  At 
this  point  the  divergence  term  in  the  above  equation  becomes  relatively  unimportant  and  can  be 
dropped.  The  resulting  system  of  ordinary  differential  equations  can  now  be  advanced  with  a 
much  larger  time  step.  This  later  regime  will  be  called  the  late  time  regime  and  is  characterized 
by  the  transition  to  an  arc. 

We  now  need  to  specify  the  various  terms  in  these  equations.  We  will  begin  with  the  mobility 
terms. 
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Electron  and  Ion  Mobility 


The  drift  velocity  of  electrons  at  equilibrium  is  found  from 


doe 

dt 


^E-v:Ue=0 

me 


resulting  in 


el-  - 

°e  = - ^E  =  PeE  • 

m  v. 


(10) 


(ID 


We  will  assume  that  the  ionization  level  is  low  enough  so  that  only  collisions  with  neutral 
molecules  need  to  be  considered  and  so 


e  1  _  e  1 
mvl  m  nQ  (crmwe) 


From  Kline  [Ref.  1]  we  find 

He(.cin 1  /  V -s)  =  2.7 x  105  / P(Torr)  . 


From  Vitello  [Ref.  7]  we  find 

jUe(cm2 /  V - s)  =  2.9x  105  / P(Torr ) 
jup(cm2  /V -s)  =  2.6xl03  / P(Torr ) 


(12) 


(13) 


(14) 


The  assumption  that  the  drift  velocity  is  proportional  to  electric  field  strength  breaks  down  for 
both  high  electric  field  strengths  and  high  ionization  densities.  Fongmire  [Ref.  30]  gives 


MePX™2  /V~s)  =  \ 


0.79 

0.25^3  x\04pr/E, 
0.079 


,Et pr  < 3x103F / m 
3x  103  <  E / pr  <  3xl05F/m  > 
,E I pr  >  3x  105F/ m 


(15) 


where  the  relative  air  density  defined  at  low  pressure  is  pr  =  no  /  N0 .  We  will  further  assume 
that  pr  is  fixed  by  the  initial  conditions,  consistent  with  the  weak  ionization  assumption  used 

here.  Vitello’s  result  for  the  electron  drift  velocity  in  nitrogen  should  also  be  modified  for  high 
field  levels.  A  possible  modification  is  (MKS  units) 
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-2 


(16) 


VePri™2  !V~  S)  = 


3.8x  10  2 
3.8x10^1. 14x10 1  prIE(VI  in) 


,EI pr  <1.14xl07F/m 
,E! pr  >1.14x107F/»z 


We  will  use  the  last  two  results  for  the  electron  mobility  in  air  and  nitrogen  respectively.  For  the 
ion  mobility  we  will  use  values  for  air  and  nitrogen  from  Longmire  [Ref.  28]  and  Vitello  i.e. 

u  p  (cm2  /V -s)  =  2.4xl0~4  -  Air 

.  (17) 

PpPr (cm  1 IV -s)  =  3 .■ 4x10  4  - Nitrogen 


Electron  Diffusion 

From  Vitello  we  get  in  MKS  units 

DePr(m2  / s)  =  0.2  .  (18) 


Electron  Collision  Ionization 


The  electron  impact  ionization  or  first  Townsend  coefficient  is  also  given  by  Vitello 


a(m  1 )  =  57()P(Torr)cxp 


r -26000P(Torr)^ 


E(V/m ) 


=  4.33x10 5 pr  exp 


-1.98x10 


7  A 


E(V  /  in)/  pt 


rj 


(19) 


We  can  now  form  the  electron  avalanche  or  cascade  rate  as  the  product  of  the  first  Townsend 
coefficient  and  the  electron  drift  velocity.  We  have 


C(sec  l^-=a[rn  1  )|t>e|(m/sec) 
=  4.33xl05/?r  exp 


r  -1.98x  107  A 


E{V  I  m)l  p r 


Me  E  ■ 


(20) 


Inelastic  Excitation  and  Decay 

The  macroscopic  cross  section  8(cm  x)  for  inelastic  excitation  of  neutral  molecules  is  given  as  a 
ratio  8 1  a  by  Yoshida  [Ref.4]  as 
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1 


s_ 

a 


Ka  J 


[1  HP /Pa)] 


Sa  A1A1  268  353.4x  102  547.3xl04 

—  =  0.101  +  - - - - ;— + - — 

a  {E/P )  {E/P)  {E/P) 


160.5  xlO6  212.4xlQ7 
{E/P)4  +  {E/P)5 


(21) 


where  P0  =  60  Torr  is  the  quenching  pressure  of  the  inelastic  transition  of  interest.  The  decay 
time  of  the  transition  is  also  given  by  Yoshida 


Tm  1 +  P/P0 


(22) 


where  ra  =  36ns  is  the  low  pressure  value.  To  be  consistent  we  will  use  from  Yoshida 


a 


(cm  1  j  =  7 ,56P(Torr)cxp 


' -212^ 

KJJpy 


(23) 


In  the  above  equations  E  /  P  has  units  of  volts/cm-Torr.  We  can  now  define  the  rate  coefficient 
for  creation  of  excited  neutral  molecules  by 


f?(sec  x^  =  S{m  1  )|r>e|(/n/sec) . 

The  rate  coefficient  for  decay  of  the  excited  molecules  is  given  by 

fi  (sec  1 )  =  — . 

V  z\„ 


(24) 


(25) 


Electron  Attachment 


For  the  electron  attachment  rate  in  air  we  will  use  the  simple  three  body  rate  from  Longmire 
[Ref.  31] 


v4(sec_1)  =  108/?2 .  (26) 

The  difference  between  attachment  and  avalanche (C-  A  )  will  determine  the  electron 
avalanche/breakdown  threshold  in  air. 

El ectron -Ion  Recomb inati on 

he  electron-ion  reaction  is  assumed  to  go  by  radiative  recombination  and  can  be  approximated 
from  Longmire  as 


Rei=  2x1(T13 


(27) 
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Ion-Ion  Recombination 


The  ion-ion  recombination  reaction  is  assumed  to  be  a  three  body  mechanism  and  can  be 
approximated  from  Longmire  as 


Ru  (/n3/s)  =  2xl  (T12  pr . 


(28) 


UV  Photo-Ionization 


In  order  to  complete  the  model  we  need  to  consider  photo-ionization  of  the  ambient  gas  and 
electron  photoemission  from  the  cathode. 

Given  the  source  S  {photons  /  cm3  -  sec)  the  ionization  rate  at  a  location  R  =  0  is  given  by 


all 

space 


r>Sp(R')^e~^dV 

4xR'2 


(29) 


or  as  given  by  Penny  [Ref.  6] 


all 

space 


hD(R')y/(R'P)PdV 
4  ttR'2 


(30) 


Penny  has  semi  log  plots  of  y/(RP )  vs  x  =  RP  for  nitrogen  and  air.  In  order  to  use  these  results 
here  the  coefficient  y/(x)  was  fit  by  a  sum  of  exponentials  of  the  form 

y/{x)  =  axe~b'x  +  a2e~hlX  + ...  (3 1) 


where  x  is  the  range  in  units  of  cm-Torr.  A  simple  fitting  procedure  was  used  where  the  plot 
was  divided  up  into  sections  and  each  section  was  fit  by  a  single  exponential  using  the  relations 


b  _  ln(i//|)-ln(^2) 

x2-xx 

ax=y/2l  e~h'x . 

For  air  we  have  a  three  term  fit  where 

ax  =  5.0x  10-4  bx  =2.705x10  1 
a2  =6.0x10  6  h2  =1.455x10  2 
a2  =  6.0 x  10”7  h3  =6.383x10  3. 


(32) 


(33) 
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For  nitrogen  we  have  a  three  tenn  fit  where 


ax  =  l.OxlO-3  ^  =3.838x10  1 

a2  =  3.0x  1CT7  b2  =  2.906x10"  (34) 

a2  =  2.0x  1(T8  =2.302x10  3. 


A  list  is  given  in  Table  1  which  gives  the  values  for  the  empirical  coefficient  y/(RP)  with  units 
of  electron-ion  pairs  created  per  ionizing  collision-str-cm-Torr  as  a  function  of  the  range- 
pressure  product.  Also  shown  is  the  effective  photon  absorption  coefficient  /i(cmlTorr _1)  for 
air  and  nitrogen. 


Table  1. 

Ionization  and  Absorption  Coefficients 


RP(cm-Torr)  y/(air )  {//(nitrogen )  /u(air)  //(nitrogen) 


1.000E+01 
4.105E+01 
7.21 1E+01 
1.032E+02 
1.342E+02 
1.653E+02 
1.963E+02 
2.274E+02 
2.584E+02 
2.895E+02 
3.205E+02 
3.516E+02 
3.826E+02 
4.137E+02 
4.447E+02 
4.758E+02 
5.068E+02 
5.379E+02 
5.689E+02 
6.000E+02 


3.919E-05 

3.771E-06 

2.480E-06 

1.648E-06 

1.106E-06 

7.508E-07 

5.162E-07 

3.600E-07 

2.550E-07 

1.835E-07 

1.341E-07 

9.963E-08 

7.510E-08 

5.738E-08 

4.439E-08 

3.470E-08 

2.737E-08 

2.176E-08 

1.741E-08 

1.400E-08 


2.178E-05 

1.093E-07 

5.385E-08 

3.074E-08 

2.076E-08 

1.613E-08 

1.373E-08 

1.226E-08 

1.120E-08 

1.034E-08 

9.590E-09 

8.914E-09 

8.293E-09 

7.719E-09 

7.185E-09 

6.689E-09 

6.228E-09 

5.798E-09 

5.398E-09 

5.026E-09 


6.331E-02 

2.463E-02 

1.903E-02 

1.623E-02 

1.446E-02 

1.319E-02 

1.224E-02 

1.147E-02 

1.085E-02 

1.032E-02 

9.872E-03 

9.481E-03 

9.136E-03 

8.829E-03 

8.554E-03 

8.305E-03 

8.079E-03 

7.872E-03 

7.681E-03 

7.504E-03 


2.239E-01 
3.331E-02 
1.558E-02 
9.61  IE-03 
6.739E-03 
5.089E-03 
4.034E-03 
3.309E-03 
2.784E-03 
2.389E-03 
2.082E-03 
1.838E-03 
1.640E-03 
1.476E-03 
1.339E-03 
1.222E-03 
1.122E-03 
1.036E-03 
9.602E-04 
8.938E-04 


UV  Photo-Emission 


In  order  to  detennine  the  electrical  discharge  UV  induced  photoemission  at  the  cathode  we  need 
the  incident  directed  photon  flux  at  the  cathode  surface.  This  flux  is  given  by 
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(35) 


rrAe-Cl)Sn(R')e^R'dV 

g-fff  1  y - 

half  ‘-tsiiv 

space 

where  en  is  the  unit  normal  to  the  cathode  surface  and  Q  is  the  element  of  solid  angle.  Again  it  is 
interesting  to  consider  the  case  where  Sp  ( R'[  is  uniform  over  the  half  space.  Rewriting  in 
spherical  coordinates  we  have 


g.  e  t=2  cQs  es  -*x R' 2dR>  sin  0d0d(/) 

Ql-  1  J  1 - £ 

i?'=0  6=0  </>=0 


4  tcR 


1 2 


o  R'= oo 

- 1  - 

■*  L 


AD’ 


dR'  = 


4/f 


The  resulting  one-sided  electron  flux  is  given  by 

T+  [el cm 2  -sec^  =  Yp  [el p)Q+p  (pi cm2  -sec). 
Explicitly  putting  in  the  pressure  tenn  results  in,  [Ref.  4]  [Yoshida,  (1),  (2)  71], 

Yr{e,-n)Sp(R')e-{'-,F>,rFdV 


ft  = 


Ilf 

half 

space 


4  7tR' 


(36) 


(37) 


(38) 


The  advantage  of  this  fonnulation  is  that  the  pressure  normalized  absorption  coefficient 
pJP  can  be  derived  from  experiment.  However  we  will  still  need  to  find  a  the  photon  source 

term  Sp  (7f )  and  the  cathode  electron  yield  coefficient  Yp  [el  p) . 


Penny  [Ref.2]  has  provided  log-log  plots  of  p(x)  =  pa!  P  for  nitrogen  and  air.  In  order  to  use 
their  results  here  the  coefficient  y/(x)  was  fit  by  a  power  sum 

ju(x)  =  alx~b'  +  apCbl  + ...  (39) 

where  //  has  units  of  cm-1  torr-1  and  x  is  the  range  in  units  of  cm-Torr.  A  simple  fitting 
procedure  was  used  where  the  plot  was  divided  up  into  sections  and  each  section  was  fit  by  a 
single  exponential  using  the  relations 


b  _  In  (a)~  In  (//2) 

ln(x2)-ln(x1)  (40) 

=  02  (x)  /  xf  . 
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For  air  we  have  a  two  term  fit  where 


ax  =  10.0  bx  =  2.736 
a2  =  0.123  b2  =0.4372. 

For  nitrogen  we  have  a  one  term  fit  where 

a,  =  5.0  hj  =1.349. 

A  list  of  values  is  also  given  in  Table  1. 

As  discussed  previously  the  photon  source  is  given  by 


(41) 


(42) 


5 


p 


(43) 


From  Yoshida  [Ref.  4]  we  find  that  the  cathode  electron  yield  is  related  to  the  first  Townsend 
coefficient  a(cm  '),  the  inelastic  excitation  coefficient  8  (cm  1 )  (macroscopic  cross  section  for 
electron  collision  induced  inelastic  excitations  of  neutral  molecule  states  that  then  decay  by 
photoemission),  and  the  second  Townsend  coefficient  y  by  the  relationship 


SYp  =  ay.  (44) 

The  second  Townsend  coefficient  y  is  a  function  of  the  cathode  material,  the  ambient  gas,  and 
the  electric  field  to  pressure  ratio  E  /  P  and  is  defined  as  the  number  of  secondary  electrons 
produced  at  the  cathode  per  ion  pair  produced  in  the  gas.  Values  typically  vary  between  1 0  and 
10  3 .  The  first  and  second  Townsend  coefficients  are  combined  to  give  the  simple  electrical 
breakdown  criteria  for  small  gaps  i.e. 


yead  =  1.  (45) 

For  electrical  discharges  in  nitrogen  at  P=300  Torr  Yoshida  also  suggests  a  value  y  ~  8.2  x  10  3 . 
The  ratio  81  a  was  discussed  previously.  Values  for  a  /  P ,  a  ,  8/  a  ,  and  8  are  given  in  Table 
2.  For  reference  5 MV  /  m  at  760  Torr  corresponds  to  66  V/cm-Torr.  Choosing 
E  /  P  =  50V  /  cm  -  Torr  results  in  8  /  a  =  1.19  and 


Yp(e!  p) 


yS  j 


r  = 


1.19 


(8.2x10  3)~  lx  10~2. 


(46) 


This  result  is  somewhat  higher  than  values  of  10  3  -10  4  anticipated  from  some  unpublished  data 
taken  by  one  of  the  authors  on  unprepared  samples  of  copper. 
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Table  2  gives  values  for  at P(cm  lTorr~l) ,  a(cml),  St  a  and  S(cm~x)  [Yoshida,  (1),  (2),  76] 
as  a  function  of  E  /  P  (V  /  cm  -  Torr )  at  a  pressure  of  760  Torr.  The  decay  time  for  the  level  is 
t  =  2.63 ns  . 


Table  2 

Electrical  Coefficients  for  Air 
-T)  ^(cmlTorrl)  a  (cm1)  St  a 


E 

—  (V / cm  - 

P 

1.000E+01 

2.000E+01 

3.000E+01 

4.000E+01 

5.000E+01 

6.000E+01 

7.000E+01 

8.000E+01 

9.000E+01 

1.000E+02 


1.163E-1 1 
9.378E-06 
8.728E-04 
8.420E-03 
3.281E-02 
8.123E-02 
1.552E-01 
2.523E-01 
3.681E-01 
4.980E-01 


8.841E-09 

7.127E-03 

6.633E-01 

6.399E+00 

2.493E+01 

6.174E+01 

1.180E+02 

1.917E+02 

2.798E+02 

3.785E+02 


7.563E+02 

1.975E+01 

4.517E+00 

2.069E+00 

1.187E+00 

7.636E-01 

5.307E-01 

3.913E-01 

3.027E-01 

2.435E-01 


8 (cm  '  ) 


6.687E-06 

1.408E-01 

2.996E+00 

1.324E+01 

2.960E+01 

4.714E+01 

6.261E+01 

7.504E+01 

8.468E+01 

9.215E+01 


Air  Heating  Model 

The  early  time  equilibrium  streamer  dynamics  assume  that  the  electron  temperature  is  directly 
related  to  the  local  electric  field  and  air  density.  In  fact,  the  molecular  gas  will  heat  up  from 
electron-neutral  collisions.  If  we  consider  the  contributions  due  to  Joule  heating  from  electron 
and  ion  currents 


J  =  e[neve  +  Yj no  )  +  \e\ £ n+u+  (47) 

and  an  additional  term  from  the  intrinsic  thermally  induced  electrical  conductivity  <jpe  of  the 
neutral  gas  we  can  calculate  the  change  of  internal  energy  from 

- igfJ  +  gE2-^  (48) 

dt  ip 

where  e  is  the  specific  internal  energy.  The  last  term  in  the  above  equation  represents  the 
radiative  loss  of  energy  from  an  optically  thin  gas  where  ahb  =  5.67  x  10  HW  /  m 1  -  deg  if4 .  The 
Plank  mean  free  path  i  p  was  not  available  so  the  Rosseland  mean  free  path  i  R  ~  10m  at  a 
temperature  of  10,000  degrees  K  was  used  Brode  [Ref.  32],  For  the  present  time  we  will  ignore 

3 

hydrodynamic  effects  and  assume  that  the  gas  density  is  constant  so  po  =  1.17688  kg/m 


18 


d_±s±  [E.J+ 

dt  p0  l  £r 


A  relatively  simple  model  for  the  electrical  conductivity  of  heated  air  has  been  proposed  by 
Plooster  [Ref.  27]  that  consists  of  a  caloric  equation  of  state  that  can  be  inverted  to  get  the 
temperature  as  a  function  of  specific  energy  density  and  an  electrical  conductivity  model. 
Plooster’ s  equation  of  state  is  written 


e(p,T)=  -(5  +  4,)  +  3(4+4)  RT  +  AJ0+Atl,+A,J2 

V  2. 


where  a  is  the  specific  internal  energy,  Ao  represents  the  fraction  of  the  air  molecules  that  are 
dissociated,  and  AX,A2  represents  the  fraction  of  air  atoms  that  have  been  signally  or  doubly 
ionized.  These  coefficients  are  functions  of  temperature  T ,  and  density  p  .  The  electron  density 
due  to  thermal  ionization  is  given  by 


nep=2PrK(Al+A)- 


Given  the  change  in  internal  specific  energy  it  is  possible  to  numerically  solve  the  equation  of 
state  for  temperature.  Then  the  electron-neutral  and  electron-ion  mobilities  are  given  by 


Men(m2/V-S)  = 


2.45xlCT 


pr{\-A)Jr{K)  p,\\-A)^v 


( m 1 IV  -s^-. 


PrA  111 


1.89x10  9r(^)3/2 

,  ( 1.63x10 -4T(K) 

L  in  - i  ;,v — - 

l  (APrf 


2.36xlQ-37;3;2 

a  i  f  L89Tev 

PrA  ln  T- - 7U3 

U  APr) 


The  electron-neutral  and  electron-ion  conductivities  are  given  by 


aenp=\e\nepPenp  =  Fl(nep’T) 
ip=\e\nePMeip=F2(neP’T) 


<J  -  = 

eip 


and  the  resulting  conductivity  is  given  by 


cr  cr  . 

en  ei 

cr  +  cr  . 

en  ei 
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Electrostatic  Model 


In  addition  to  creating  and  transporting  charge  it  is  necessary  to  advance  the  electric  fields  in  a 
self-consistent  manner.  The  approach  used  here  is  to  advance  Poisson’s  equation 


dr 2  r  dr  dz  so 


(56) 


where  the  net  charge  density  at  the  corners  of  the  finite  difference  cell  are  detennined  by 
summing  the  free  electrons,  positive  ions  and  negative  ions  i.e. 

P  =  e(ne  +  YJn-)  +  NIX-  (57) 


The  electric  field  is  then  found  from  E  =  -V <f> .  Solving  Poisson’s  equation  also  requires  the 
electron  boundary  conditions.  For  the  early  time  discharge  we  will  define  an  effective  discharge 
current 

zH  r  r  \  z„ 

Id  =  I  |  J.  7.  nr  dr  dz  1 J  dz.  (58) 

Z,  V  0  J  z, 

The  voltage  drop  across  the  electrodes  is  then  given  by  Vu  -  V,  =  Vs  -  IdRs 
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3.0  Numerical  Formulation 

The  coordinate  system  is  shown  in  Figure  5.  The  outer  boundaries  are  defined  on  the  lower, 
upper,  and  outer  surfaces  by  boundary  potentials  VL,VU,VO . 

</>(i,l)  =  VL 

<f>(i,nez)  =  VU  (59) 

<j)(ner,  k )  =  VO. 


The  object  potentials  Vt  (pot_modmap(i,k))  are  also  set  along  with  any  fixed  charge 
P fixed  tt’k)  and  the  volume  charge  p(i,k)  is  set  to  zero.  A  direct  solution  technique  described  by 
Kunhardt  [Ref.  25]  ,  and  Hockney  [Ref.  26]  was  adapted.  We  begin  by  with  Poisson’s  equation 

P_ 

£„ 


dr"  r  dr  dz 


(60) 


and  expand  the  potential  and  charge  density  as  discrete  sine  transforms  in  the  z  direction 

oo 

<!>(r,z)  =  YJ<j>n(r)sm 


f  \ 
nnz 


n= 1 


P(r,z)  =  Yjpn(r)  sin 


f  \ 

nnz 


(62) 


V  ^ma x  J 


Note  that  this  expansion  assumes  that  the  potential  is  zero  at  the  lower  and  upper  outer 
boundaries.  This  is  not  really  a  limitation  as  we  can  place  inner  surfaces  of  arbitrary  potential  to 
build  the  desired  problem.  Using  this  approach  along  with  the  assumption  of  a  uniform  mesh  in 
the  axial  or  z  directions  it  can  be  shown  that 


Figure  5.  Numerical  grid  for  the  two  dimensional  problem. 
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(63) 


aV„(r)  |  1  a^,(r)  ,  2 

d2r  r  dr  Az2 


(  \ 

nn 

cos 

U+iJ 

(64) 


This  equation  is  tridiagonal  on  the  2D  finite  difference  mesh.  The  method  of  solution  is  as 
follows.  First  the  quantities  p„(r)  are  detennined  by  fast  sine  transforms  as  described  by  Press  et 

al  [Ref.  33].  Then  we  solve  the  tridiagonal  equation  shown  above  for  (f>n{r)  •  Finally  we  find  the 
potential  on  the  two  dimensional  mesh  by  the  fast  sine  transform. 


=  sin 

n= 1 


(  \ 
nnz 

V  ^ m ax  / 


(65) 


In  order  to  use  this  direct  approach  it  is  necessary  to  convert  the  electrical  potential  of  defined 
objects  or  surfaces  within  the  grid  to  equivalent  electrical  charges.  First  consider  an  interior 
potential  surface  as  shown  in  Figure  6. 


Figure  6.  Placement  of  a  potential  surface  in  the  finite  difference  grid. 
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We  begin  by  placing  a  unit  charge  (ql  =  1  /  e0 )  at  the  first  object  node  on  the  potential  surface  and 

solve  Poisson’s  equation  using  the  direct  solver  to  find  the  potentials  at  all  the  nodes  on  the 
potential  surface.  We  then  set 


au  Vx 

ai\  ~  ^2 


(66) 


a.i=V.. 


We  then  repeat  this  process  for  the  rest  of  the  nodes  and  form  a  system  of  equations  which 
relates  the  potentials  at  the  object  nodes  to  the  charges  at  the  nodes  i.e. 


or 


Vi\ 

ai  i 

a\  2 

•  <hn 

9, 

V2 

= 

a2\ 

a2  2 

■  a2n 

q2 

(67) 

k\ 

G,2 

■  ann_ 

3n'_ 

V  =  C  xq.  (68) 


As  a  result  by  calling  the  direct  solver  n  times  we  can  form  the  inverse  capacitance  matrix  C  1 . 
This  matrix  is  normally  not  excessively  large  and  can  be  inverted  to  get  the  capacitance  matrix 

C .  This  operation  only  has  to  be  done  at  the  beginning  of  the  calculation  and  only  if  the 
geometry  changes.  The  spark  code  checks  to  see  if  the  capacitance  matrix  needs  to  be 
recalculated  at  the  beginning  of  each  run. 

Given  the  capacitance  matrix  the  Poisson  calculation  proceeds  as  follows.  First  the  direct  solver 
solves  for  the  potentials  given  the  ionization  charge  but  ignoring  the  boundaries.  The  resulting 
(false)  potentials  V*  at  the  object  nodes  are  then  used  to  calculate  false  charges  via 

q*=dv*.  (69) 

Then  the  negative  of  the  false  charges  is  added  to  the  ionization  charge  and  the  direct  solver 
again  solves  for  the  potentials.  At  this  point  all  the  potentials  at  the  object  surface  are  zero.  If 
the  object  potentials  are  not  zero  we  only  need  to  subtract  them  from  V*  and  solve 

r=S(v-r^y  (70) 

Once  the  potentials  are  determined  the  initial  electric  fields  are  found  from 
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er(i,  k)  =  -  [0(z  +  1, A:)  -  OO,  A:)]  /  Ar(/) 
ez(z,  A:)  =  k  + 1)  -  0(z,  A:)]  /  Ax(k). 


(71) 


We  also  need  to  initialize  the  late  time  Plooster  equation  of  state  variables.  The  only  difference 
is  adding  an  additional  dimension  to  the  previous  results.  Given  the  initial  neutral  gas 
temperature  T (z,  k)  and  density  po  we  have 


se{i,k)  =  s{p0,T{i,k )) 

(72) 

(/, k)  =  2 prNa  ( 4 (T 0,  A :))  +  4 (T (z,  A:))) 

(73) 

(h  *)  =  (h  *)> T Xh  kj) 

°enP  (Uk)  =  F2  (z nep(i,k),T(i,k) ) 

(74) 

M  venp(kk)(j  (i,k) 
cr  (i,k)  = 

venp(i,k)  +  aeip(i,k) 

(75) 

For  the  nth  time  step  advance  from  timet"  to  time  t"  +  At  =  t"+l  we  begin  with  the  number 
densities  n"  (z,  k),  n"+  (/,  k),  n"  (/,  k),  n"(i,  k)  for  the  ionization  electrons,  positive  ions,  negative  ions, 
and  neutral  molecules.  We  also  have  the  electric  field  er"  (i,k),ezn  (i,k) .  We  will  initialize  the 
final  electric  field  ern+1  ( i ,  k),  ez""  (/,  k )  to  er"  ( z ,  k),  ez"  ( i ,  A:)  and  then  iterate  until  the  change  in 
er"+1  (z,  A:),  ezH+1  (z,  A:)  is  acceptable.  Therefore  we  start  the  iteration  by  setting 

ern+l(i,k)  =  er"(i,  k ) 
ezn+1(i,k )  =  ez"(i,k) 

(76) 

ez'  _  temp(i,  k)  =  er"  (z,  A:) 
ez  _  temp(i,  k)  =  ez"  (z,  k). 

We  now  loop  over  the  j  index.  The  first  step  is  to  define  time  averaged  fields 


er(i,k)  =  0.5^ern(i,k)  +  ern+l(i,k)^j 
ez(i,  k )  =  0.5  (ez"  (z,  k )  +  ez"+1  (z,  A')) . 


(77) 


These  fields  are  not  located  at  the  fluid  cell  center  but  rather  at  the  edges.  Therefore  we  take 
spatial  averages 
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(78) 


( er(i,k ))  =  0.5(er(i,k)  +  er(i-l,k)) 
(ez(i,  k)}  =  0.5  (ez(z,  k )  +  ez(i,  k  - 1)) 


emag(i,k )  =  (er(i,k)f  +  (ez(i,k)y 


at  the  fluid  cell  centers  ( x(i),z(k )  )  and  calculate  the  mobilities,  rate  coefficients,  and  drift 
velocities  as  was  done  for  the  one  dimensional  case  and  drift  velocities  as  was  done  for  the  one 
dimensional  case.  For  air  we  have 


MUk) 


1 


Pr 


0.79  ,  emag(i,  k)/  pr  <3x10 3  F  /  m 

<  0.25^3  x  10 4  pr  /  emag(i,  k )  ,  3  x  1 03  <  emag(i,  k)  /  pr  <  3  x  1 05  V  /  m  > 

0.079  ,  emag(i,  k)  /  pr  >  3  x  1 05  V  /  m 


(79) 


pp(i,k)  =  2Ax\0  4  / pr 


(80) 


f 

a(i,k )  =  4.33xl05yOr  exp 

V 


-1.98xl07  ^ 
emag(i,  k  )  /  pr  y 


(81) 


A(i,k)  =  108/?(2 

Rei(i,k)  =  2xl0~13 

RjiiUk)  =  2x\0~12  pr 

vre  (z,  k )  =  -//e  ( emag{i ,  A:))  •  (ef(z,  A:)} 
A:)  =  -pe(etnag(i,  k ))  ■  (ez(i,  k)) 
vr+  (; i ,  A:)  =  //p  (emag(i,  k ))  •  (er.  A:)) 
yz+  (/,  A:)  =  (emag(i,  k ))  •  (ez,  A:)) 

(z,  A:)  =  -//p ( emag{i ,  A:))  •  (er (z,  A:)) 
l>z  (z ,  A:)  =  —p  ( emag(i ,  A:))  •  (ez(z,  A:)) 


(82) 

(83) 

(84) 


(85) 


C(i,k)  =  a(i,k)\oe(i,k)\ 

B{i,k)  =  8i\ve(i,k)\  (86) 

E(i,k)  =  \/rm. 
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jup(i,k)  =  3.4x10 -4  /  pr 

(89) 

a{i,k)  =  4.33x10''/?  exp  - 1-98x10 — 

{  emag(i,k)  /  pr  y 

(90) 

o 

II 

S' 

(91) 

Rei(i,k)  =  2xl0”13 

(92) 

Rii(i,k)  =  2x\0-Upr. 

(93) 

As  in  the  one  dimensional  case  we  need  to  determine  the  UV  electron  sources  and  the  UV 
induced  photoemission  at  the  cathode.  For  the  UV  electron  sources  we  will  assume  that  only 
electron  collisions  on  the  axis  are  important.  As  a  result  we  can  write 

kmax 

Se  ( i,k )  =  A  £  (hD(l,k'))Smt  ( i,k,k ') 

k'= 1 

(94) 

where  Se  ( i,k )  has  units  of  electrons  per  unit  volume  per  unit  time  and 
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,_,v  W(R'P)-P. AV(l,k') 

int '  ’  ’  '  4;rA'2 

74 2  =  r(i)2  +  (z(k)  -  zoiki))2  (95) 

AV(l  ,k')  =  7zro(\)2  Az(k') 

(nD(l,k)}  =  0.5(C(l,k)ne(l,k)  +  C(l,k  +  \)ne(l,k  +  l)). 

The  A  coefficient  is  an  empirical  scale  factor  to  account  for  the  fact  that  the  discharge  has  a 
radius  larger  than  one  cell.  It  would  be  expected  that  A  ~  the  number  of  radial  cells  in  the 
discharge. 

For  the  UV  induced  photoemission  we  first  determine  an  electron  emission  flux  at  the  cathode 
defined  by 


/ctnax 

K(hA)  =  Y^e/p)B2(e,-n){Sp(l,n)ri,,(i,^U')  W 


k'= 1 


where  S  (l  ,k')  has  units  of  photons  per  unit  volume  per  unit  time  and 


YmX(ic,kc,\,k')  = 

(v^)  = 


-M(R'P)R'P 


AV(l  ,k') 


AkR'2 

\zo(k')-z(kc) 


Wc) 


1 


,4  >1 


R'2  =r(ic)  +  (z(kc)-zo(k')) 
AV(l,k')  =  7rro(l)2Az(k') 


Sp(lk')  = 


(97) 


The  first  term  in  the  sum  represents  the  en  ■  Q  cosine  tenn  in  the  integration  over  the  emitting 
surface.  The  units  of  the  absorption  coefficient  are  length  1  pressure  1 .  As  a  result  the  integral 
Yint  (/),/')  has  units  of  length.  Given  the  electron  emission  flux  the  electron  source  at  the  cathode 
is  given  by 


S.(ie,ke)  = 


r:a.*c) 


zo(kc)  -  zo(kc  - 1) 


-.th 


(98) 


As  in  the  one  dimensional  case  after  the  air  chemistry  and  UV  sources  are  called  the  electron  and 
ion  transport  tenns  are  detennined  using  a  second  order  upwind  differencing  scheme.  We  define 
the  advective  fluxes  as 
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Yxadv(i,k)  = 

TZaAUk)  = 

We  also  include  the  diffusion  terms 


\ n(i,  k )  (vr(i,  k )) ,  ( ur(i ,  k ))  >  0 
n(i  + 1,  k )  (n>r(z,  A;)) ,  (vr(i,  A:))  <  Oj 
f zz(z,  A:)  (t>z(z,  A:)) ,  (uz(z,  A:))  >  0 
[zz(/  + 1,  A:)  (n>z-(z,  A:)) ,  (n>r(z,  A:))  <  Oj 


r rdiff  (*>  *)  =  - ( «" O'  +  hk)  -  n"  (z,  A:)) 


Fz/,y/  (0  *)  =  ~  - “77  (""  (*>  k  +  !)  -  n”  (i,  A:)) . 


Wv’  '  A z(k) 

The  transport  terms  are  now  advanced  using 


(99) 


(100) 


_v-(zzu)]+[v-rrf#_  = 


1  \^racIv{i,k)AR(i,k)-Tradv(i-\,k)AR(i-\,k),  vr(i,k)>  Ol 


(101) 


A V (z, k )  [r radv (z  + 1, k)AR(i, k)-T radv (z, k)AR(i -\,k),  ur(i, k )  <  Oj 

_ I _ [Tzadv{i,k)AZ(i,k)-Tzadv(i,k-\)AZ(i,k-\),  uz(i,k)>()\ 

A V(i, k )  [rza(/i,(z, k  +  l)dA?(z, k) ~Tzadv(i, k)AR(i, k - 1),  vz(i, k )  <  0 J 

+  AF|.  —  fF^,//  0,  k)AR(i,  k )  -  ix^  (z  - 1,  k)AR(i  - 1,  A:)) 

+  [Tzdiff  (O  ^)ZZ0>  *)  -  Fztz«r  (h  k  ~  ])AZ(k  k-\j). 

Given  the  advective  flux  the  current  density  is  given  by 

cur(i,k)  =  eTradv(i,k) 
cuz(i,k)  =  eYzadv(i,k). 

the  change  in  number  densities  due  to  advection  and  diffusion  must  be  combined  with  the  growth 
or  decay  terms  due  to  avalanche,  attachment,  recombination,  and  so  on.  Once  the  number 
densities  are  advanced  the  early  time  electrical  conductive  is  calculated  from 


(102) 


<t(z,  k)  =  \e\  (/4(z,  k)nf  (z,  k)  +  pp  (z,  k)  ( nn++l  (z,  k )  +  z?"+1  (z,  A:))). 


(103) 


This  conductivity  is  not  directly  used  in  the  calculation  but  serves  as  a  diagnostic.  Next  the 
electrode  and  wall  boundary  conditions  are  set  using  the  modmap  (i,k)  array.  The  net  charge 
density  is  detennined  from 


pn+l  (z,  k )  =  e[f z"+1  (z,  k )  -  <+1  (z,  k )  +  zz"+1  (z,  A:)) 


(104) 


and  the  total  net  charge  is 
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(105) 


a", =I2>"'(U)  m,k\ 

k= 1  i=l 

In  order  to  account  for  the  resistance  of  the  power  supply  it  is  necessary  to  calculate  an  effective 
discharge  current  between  the  electrodes,  assumed  to  be  located  at  grid  points  id  _  lower  and 
id  _  upper  respectively.  We  have 


kd  upper- 1  ner 

X  I  cuz(i ,  k)AZ  (i ,  k)Az(k ) 

kd  _  lower  i= 1 

kd  upper -1 

kd  lower 


(106) 


Now  we  will  adjust  the  potential  at  the  electrons  using  the  source  resistance  of  the  power 
supply  i.e. 


<f>(id  _upper)-(j)(id  _lower)  =  (Vidupper-VidJower)-RsId. 


The  Poisson  equation  solver  described  previously  is  called  to  solve 


VI  in+\  [  ii+l  n+ 1  .  ii+l  t 

9  =  — ( ne  -n+  +  n_  ) 

£„ 


(107) 


(108) 


and  the  fields  are  advanced  as 


—n+ 1  V7  J.n+1 

e  =  -V  (p  . 


We  now  check  the  convergence  of  the  iteration  by  checking 

\i,k)-e  _temp{i) | 


max 


en+l‘ 


en+l  (0  +  1-0x10 


-12 


<0.001 


(109) 


(110) 


or  noting  that  the  number  of  iterations  j  >  ymax  .  If  the  exit  criteria  are  not  met  we  set 

er  _  tempi},  k)  =  er"'x  ( i ,  k ) 
ez  _  temp(i,  k )  =  ez"+1  (/,  k) 


(HI) 


and  repeat.  If  the  exit  criteria  are  met  we  proceed  and  detennine  the  number  of  electrons  and 
positive  ions  per  unit  area  i.e. 
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(112) 


^;+1=II<+1(a)Anu) 

1  1 

nez  ner 

K+1  ^^nfii^AV^k). 
1  1 


This  completes  the  early  time  advance  where  propagations  effects  are  important  and  where  the 
electron  temperature  is  not  in  equilibrium  with  the  neutral  molecules.  However,  we  need  to 
advance  Plooster’s  equations  for  the  neutral  gas  specific  energy,  temperature,  and  conductivity  in 
order  to  know  when  to  switch  over  to  the  late  time  model.  The  specific  energy  is  advanced  by 
considering  the  thermal  energy  density  with  contributions  to  both  the  early  time  electron  flow 
and  the  heating  from  the  Plooster  conductivity.  Specifically  we  advance 

(113) 

-E)  =  0.5  ((er(/  - 1,  k)cur(i  - 1,  k)  +  er(i,  k)cur(i,  k))  +  (ez(i,  k  -  X)cuz(i,  k- 1)  +  ez(i,  k)cuz(i,  k)[) 

( ER_ave )  =  0.5(er(i-\,k)  +  er(i,  k)) 

[EZ  _  ave )  =  0.5  (ez(i,  k  - 1)  +  ez(i,  k)) 

( ps)  =  posn  (/,  k)  +  ( J  •  is)  At  +  <jep  (/,  k)[(^ER  _  ave)"  +  [EZ  _  avef  j  At 
£n+\i,k)  =  (ps)l  pa. 


In  order  to  detennine  the  temperature  we  need  to  invert  Plooster’s  equation  of  state.  The  method 
for  determining  temperature  goes  as  follows.  Given  the  specific  internal  energy  density  s" "  the 
temperature  TnA  can  be  detennined  by  fonning  the  function 

g[Tn+l)  =  £(Tn+l,p)-£n+l  =0.  (114) 


Then  by  Newton-Raphson  iteration 


r+1(y  +  i)  =  r+1(y)- 
=  T"+1{j) 


g(T‘:t'U)) 

dg(r*'(j))/dT 

(■ s{rM(j),p)-s"")&T 
-  e(r"(j)  +  AT,p)-e(r"(j),p)' 


(115) 


The  AT  term  can  be  set  to  about  Tn  ■  10  3 . 


The  Plooster  conductivity  aep(i,k)  is  advanced  as  described  previously.  The  resistance  array 

between  the  fluid  cell,  shown  in  Figure  7,  centers  is  determined  as  follows.  For  the  vertical 
resistors  on  the  axis  we  have  (i  =  1  ,k  =  1, nhz) 


RZ(\,k) 


zo{k)-z{k)  z(k  +  1)-  zo{k) 
nro( l)2  <jep  (1,  k )  xro(\f  <rep  (1,  k  + 1) 


(116) 


30 


^(1)  R(2)  R(3) 


Figure  7.  Representation  of  the  finite  difference  conductivity/resistivity  mesh  used  to  solve  for 
the  late  time  discharge.  Note  that  objects  with  fixed  potentials  are  represented  by  current 
sources. 

For  the  inner  vertical  resistors  (/  =  2 ,...,nhr,k  =  l,nhz )  we  have 


RZ{i  k,  _ zo(k)-z(k ) _ + _ z(k  +  \)-zo(k) _ 

7r(ro(J.f  -ro(i-l)2)crep(i,k)  n(ro{if  -ro{i -\)2)<Jep{i,k +  \) 

For  the  outer  vertical  resistors  (/  =  ner,k  =  1 ,  nhz)  we  have 


(117) 


RZ(ner,  k )  =  - 


zo{k )  -  z{k) 


-  + 


z(k  + 1)  -  zo(k ) 


n 


(r(ner)2  - ro(nhr)1 ) <jep(nei\k)  7i[r(nerf  - ro(nhr)2 ) aip(nei\ k  +  \) 


(118) 


Normally  the  outer  boundary  is  conductive  so  these  resistors  are  set  to  a  small  value.  For  the 
interior  radial  resistors  (/  =  1  ,...,nhr,k  =  2,...,nhz)we  have 
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ro(i)-r(i ) 


r(i  +  1  )-ro(i) 


(  ’  ^  2nro{i)lSz{\)aep  (i,  1)  2/rro(i)Az(l)crep  (i  + 1, 1) 


RR(i,k)  = 


ro(i)-r(i ) 


-  +  - 


r(i  +  \)  -  ro(i) 


RR(i,nez )  = 


2nro(i)  ( zo(k)  -  zo{k  - 1))  oep  (i,  k )  2nr<){i){zo{k)  -  zo{k  - 1))  oep  (i  + 1,  k) 

ro(i )  -  r(i )  r(i  + 1)  -  ro{i) 


(119) 


-  +  - 


27iro(i)Az( nhz ) o  v  (i,  nez)  2^ro(i)Az(nhz)cr  (i  + 1,  nez) 


The  radial  resistors  across  the  bottom  and  top  of  the  mesh  form  the  conductive  boundary 
condition  so  they  are  also  set  to  a  small  value  as  are  the  resistances  inside  any  conductive  object 
within  the  finite  difference  mesh. 


At  this  point  we  need  to  calculate  the  voltages  at  the  cell  nodes  and  the  electric  fields.  There  are 
two  approaches  to  consider.  The  first  is  to  approximate  the  resistance  between  two  objects  or 
potential  surfaces  within  the  finite  difference  mesh  by  assuming  axial  current  flow.  In  this  case 
we  can  divide  the  problem  into  a  set  of  conductive  disks  in  series.  The  admittance  of  each  disk  is 
given  by 


i  _  upper 

Y(k)=Yj  1.0/  RZ(i,k)  (120) 

1=1 

and  the  Plooster  resistance  between  the  two  electrodes  is  given  by 

kd  upper 

Z  i.o/m).  (no 

k=kd  lower 
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4.0  Test  Simulations 


In  order  to  test  the  Poisson  solver,  a  couple  of  simple  problems  have  been  solved.  In  both  cases, 
the  program  was  applied  to  a  problem  with  128  r  grid  cells  and  256  z  grid  cells.  Both  the  r 
and  the  z  are  uniform,  with  a  cell  size  of  one  half  unit.  The  first  test  case  is  for  a  sphere  of 
uniform  charge  density  of  magnitude  1/10.  The  sphere  has  a  radius  of  five  units.  Figure  8  shows 
the  charge  density  used.  The  potential  computed  from  this  charge  distribution  is  shown  in  Figure 
9. 


Cl 
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Figure  8:  Test  Run  Charge  Density. 


Figure  9:  Test  Run  Scalar  Potential. 
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Scalar  Potential  from  Past  Poisson  Solver 


This  is  a  particularly  nice  test  case  because  the  solution  of  the  resulting  equation  is  known 
analytically.  In  particular,  consider  solving  Poisson’s  equation  solved  for  a  uniformly  charged 
sphere  at  the  origin 


vV  = 


0  <r  <a 
a  <  r. 


(122) 


In  spherical-polar  coordinates,  symmetry  requires 


d(j) 

~80 


in  the  radial  variable  only: 


d(j) 

d<p 


=  0,  and  we  are  left  with  an  equation 


j_sf 2 Ml- \p°  °-r-a 

r2  dr  \  dr  J  [  0  a  <  r. 
Multiplying  by  r2  and  integrating  with  respect  to  r  once  gives 


(123) 
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(124) 


r 


2 


d(j) 

dr 


—  r3  +  c,  0  <r<a 
<  3 

c2  a  <  r. 


Dividing  by  r  and  integrating  with  respect  to  r  again  gives 


—  r2  -  — +  c3  0 <r<a 

6  r 

c2 

— -  +  c4  a  <r. 

r 


(125) 


The  first  constant  of  integration  may  be  determined  by  a  regularity  condition.  The  potential  must 
be  finite  as  r  — >  0  which  forces  c,  to  be  zero.  We  also  require  that  the  potential  must  go  to  zero 

as  r  becomes  infinite.  This  forces  c4  to  be  zero.  A  requirement  that  the  potential  and  its  first 
(radial)  derivative  are  continuous  everywhere,  but  at  r  -  a  in  particular,  gives  two  equations  that 
may  be  solved  for  c2  and  c3. 


1  2  Co 

-p0a  +c3= — - 
6  a 


and 


]_ 

3 


p0a  = 


result  in 


c2 


2  Poa  ’  C3 


finally  giving  the  solution 


(126) 


(127) 


(128) 


4>  = 


A)r2_I 
6  2 
1  p0a 3 


p0ci2 


3  r 


0  <r  <a 

a  <  r. 


(129) 


The  scalar  potential  given  by  this  analytical  solution  for  a  sphere  analogous  to  the  one  used  in 
Figure  9  is  shown  in  Figure  10.  The  difference  between  the  analytic  solution  and  the  one 
obtained  from  the  fast  Poisson  solver  is  shown  in  Figure  11. 
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There  are  two  primary  causes  for  the  difference  in  the  values  of  the  analytic  potential  and  the 
computed  one.  The  first  is  that,  due  to  discretization,  the  charged  objects  are  slightly  different. 
This  is  clearly  seen  in  Figure  1 1  where  the  error  has  large  peaks  near  the  boundary  of  the  charged 
sphere.  The  second  source  of  error  is  due  to  the  fact  that,  for  the  numerical  solution,  the  scalar 
potential  is  forced  to  be  zero  at  a  cylinder  defined  by  (polar)  radius  r  =  R  and  with  the  top  and 
bottom  defined  by  z  =  0  and  z  =  Z,  respectively.  The  analytic  solution  asymptotically 
approaches  zero  as  the  (spherical)  radius  becomes  large,  but  is  not  zero  for  any  finite  value. 


Scalar  Potential  from  Analytic  Solution 


Difference  in  Scalar  Potential:  Analytic  - ! 


Figure  10:  Analytic  Scalar  Potential  for  a  Uniformly 

Figure  1 1 :  Difference  between  Analytic  and  Fast 

Charge  Snhere. 

Poisson  Solutions. 

To  give  the  algorithm  a  bit  more  of  a  challenge,  a  second  test  case  was  run.  In  this  case,  two 
spheres  of  equal  but  opposite  charges  are  displaced  from  the  center  of  the  grid  by  three  units, 
each  in  a  different  direction.  This  results  in  two  equal,  but  opposite  crescent  shaped  charge 
objects.  In  the  area  near  the  origin  where  the  two  spheres  overlap  results  in  a  cancellation  of 
charge.  Figure  12  shows  the  resulting  charge  density  for  this  run.  Figure  13  shows  the  resulting 
scalar  potential  computed  from  the  charge  distribution  of  Figure  12. 


Figure  12:  Charge  Distribution  for  Second  Test  Case. 


Figurel3:  Scalar  Potential  for  Second  Test  Case. 
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5.0  Conclusions 


Numerical  models  have  been  developed  for  the  range  of  physical  process  on  spatial  and  temporal 
scales  required  for  the  analysis  of  the  USP  Laser  initiated  high  voltage  discharge  problem.  These 
specifically  include: 

Treatment  of  the  reaction  rate  equtions  for  determining  the  evolution  of  the  numerous 
charged  particle  species  in  air  under  high  electric  field  strengths.  This  treatment  is  used  in 
the  early  time  discharge  initiation  phase  of  the  process. 

A  very  fast  method  for  determing  the  quasi-static  values  of  the  electic  field  through  out 
the  simulation  volume,  at  each  time  step  of  the  simulation.  This  was  accomplished  by 
development  of  an  FFT/Cubic  Spline  based  Poisson  algorithm. 

A  simulation  of  the  production  and  absorption  of  ionizing  uv  radiation  within  the 
computational  volume. 

A  method  for  representing  various  metalic  object  boundry  conditions  within  the  overall 
cylidyrically  symmetric  simulation  volume  and  still  allow  for  the  zero  field  end 
boundaries  required  by  the  FFT  base  poison  algorithm.  This  method  employs  the 
placement  of  a  fictitous  charge  distribution  that  yields  the  proper  boundry  conditions  for 
imbedded  conducting  surfaces. 

A  simulation  of  the  motion  of  electrons  within  the  simulation  volume,  this  was  treated 
with  a  two  step  fluid  description  of  the  moving  electrons. 

In  order  to  extend  the  computational  range  within  feasible  computer  simulation  time 
approximations  were  developed  for  use  in  the  early  and  late  time  phases  of  the  simulation: 

An  approximation  was  developed,  for  use  in  late  time  after  charge  flow  equilibration,  to 
treat  the  electron  motion  as  a  current  with  a  conductivity  based  on  an  caloric  eqution  of 
state  previously  developed  by  Plooster. 

Definition  of  a  procedure  to  treat  long  range  (Meter  scale)  discharges  as  a  series  of 
detailed  local  ‘snap  shot’  calculations  at  points  along  the  extended  streamer  path. 


The  above  outlined  models  were  integrated  and  demonstrated  on  a  single  CPU  PC  based 
workstation.  Test  simulations  were  done  to  assure  the  proper  operation  of  the  simulation 
modules.  Methods  for  parrallelizing  the  overall  simulation  for  porting  to  a  multi-core  High 
Performance  Computer  were  studied  and  outlined. 

The  modules  outlined  above  were  however  integrated  and  a  series  of  small  scale  test  problems 
completed.  These  calculations  showed  the  apparent  proper  behavior  of  the  simulation,  the 
expected  reduction  of  the  air  breakdown  voltage  associated  with  the  discharge,  and  the  formation 
and  initial  propagation  of  a  discharge  channel  were  observed.  Larger  scale  calculations  were  not 
feasible  on  a  workstation  type  computer.  Methods  for  moving  to  parallel  operation  on  an  HPC 
were  investigated  and  outlined. 
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